Exact exchange-correlation potential of a ionic Hubbard model with a free surface. 
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We use Lanczos exact diagonalization to compute the exact exchange correlation potential (y xc ) 
of a Hubbard chain with large binding energy ("the bulk") followed by a chain with zero binding 
energy ("the vacuum"). Several results of density functional theory in the continuum (sometimes 
controversial) are verified in the lattice. In particular we show explicitly that the fundamental gap 
is given by the gap in the Kohn-Sham spectrum plus A xc , the jump on v xc in the bulk when a 
particle is added. The presence of a staggered potential and a nearest-neighbor interaction V allows 
to simulate a ionic solid. We show that in the small hopping amplitude limit A xc = V in the ionic 
regime, while in the Mott regime A xc is determined by the Hubbard U interaction. In addition we 
show that correlations generates a new potential barrier in v xc at the surface. 



Density functional theory (Tj-Q (DFT) plays a major 
role in our understanding of ground state properties of 
materials. However most approximate DFT approaches 
fail to predict the fundamental gap Ac of insulators and 
semiconductors (band gapproblem) 04l3|, in systems 
ranging from bulk Silicon [8[ to ZnO [12| and other cor- 
related insulators [13|. 

At the heart of almost all practical computational 
schemes based on density functional theory [l|-[3j (DFT) 
lies the assumption, first introduced by Kohn and 
Sham [2], that the ground- state density p of an inter- 
acting electron gas in an external potential can be re- 
produced by a system of non-interacting electrons in an 
effective potential ^ks- The effective potential can be ex- 
pressed as the sum of three contributions: the external 
potential, the Hartree potential vh, and a term which 
accounts for exchange and correlation effects, v xc . The 
latter is the functional derivative of a universal "divine 
functional" [l4| of the density whose precise form is not 
known. As first discussed by Perdew et aJ. @, [H[ and by 
Sham and Schiilter [6|, the exchange-correlation poten- 
tial v xc may have a jump of order one when one particle 
is added to a solid. This jump, which is absent in local 
and semi- local approximate functionals,[9|, [l2| may ac- 
count for the error on the fundamental gap according to 
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FIG. 1: Structure of the system consisting of Lb bulk and 
Ly vacuum sites. The bulk sites have a large binding energy, 

W . 



the highest occupied Kohn-Sham eigenvalue and whose 
validity has also been subject of controversies^ [l8l - [22| . 
This debate along with the need to understand and cor- 
rect the deficiencies of approximate DFT approaches has 
recently revived the interest in small systems (zero di- 
mensional) whose exchange-correlation potential can be 



(i) 



where denotes the single-particle gap in the Kohn- 
Sham non-interacting system. The size of this effect has 
been however long debated 

In a pioneering work Gunnarsson and Sch6nhammer[7] 
studied a model of a one dimensional spinless insula- 
tor and found that A xc is small in the band insulating 
regime. Other authors have, however, argued that the 
discontinuity should be large and it should account for a 
large part of the band gap problem (3, [To|, [Hj. The elu- 
siveness of A xc is such that even its existence has been 
recently questioned [16]. 

Eq. ([]]) is based on the DFT version of Koopmans 
theorem^, [l8| which identifies the ionization energy with 



26| or lattice 



calculated exactly or very accurately [22 
systems where DFT or approximate lattice DFT schemes 
can be tested and analyzed in a controlled environment 
retaining many of the subtleties of the many-body prob- 
lem in extended systems which can not be accessed 
otherwise [2 7H35l | . 

In this work we use lattice DFT to investigate the 
band-gap problem. We calculate numerically the exact 
exchange-correlation potential of a correlated insulator 
described by a generalized Hubbard model which can 
be tuned continuously from a ionic to a Mott insulating 
regime [36]. We consider an open system with a free sur- 
face which removes any possible ambiguity related to the 
validity of Koopmans theorem and or Eq.[TJ i.e. we com- 
pute each term on the left and right of Eq. [T] separately 
which serves as a numerical test of the equation itself. We 
find that the contribution of the exchange-correlation po- 
tential discontinuity to the charge gap is non-negligible 
in both regimes. The presence of the surface also allows 
us to highlight the appearance of an anomaly in the exact 
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exchange correlation potential in the vacuum sites which 
appears as the system enters the Mott phase. 

We consider a Hubbard chain of Lb sites with a large 
binding energy called "the bulk" followed by a chain of 
Ly sites with zero binding energy termed "the vacuum" 
with open boundary conditions as shown in Fig. HJ The 
bulk is thus a truly open system which is crucial to com- 
pletely determine the exchange-correlation potential. 

The total Hamiltonian can be written as H = T + 
Hjj + H v with 

T = -t^2(cl a c x + la -n X(T + H.c.) 

XG 

Hu = U^n xt n xi + V^2n xa n x +i a (2) 

x xa 

Hy — ^ ^ ^x^xa •> 
xa 

where c X(J creates an electron with spin a =t, 4- at site x, 
U and V are respectively the Hubbard interaction and 
nearest-neighbor interaction, t is the nearest-neighbor 
hopping and we set n xa = c xa c xa . We included a con- 
stant energy shift in the lattice kinetic energy T so that 
single particle energies are measured from the bottom of 
the band. In order to simulate the work function of a solid 
the potential in the bulk is taken as v x = — wo + S(— l) x 
where wq is a large positive constant such that all parti- 
cles in the system are bound in the bulk region and the 
second term is a site dependent potential. The potential 
in the vacuum is by definition v x = 0. 

We apply DFT to this problem by considering the site 
occupancy p x = ^2 a (n xa ) as the fundamental variable Q. 
The charge density and the ground state energy are ob- 
tained using Lanczos exact diagonalization[37j. The ex- 
change correlation potential is obtained form the exact 
density inverting the Kohn-Sham problem [38|. 

In order to illustrate our capability obtain the abso- 
lute value of the exchange correlation potential we first 
consider the case of a constant external potential in the 
bulk (5 = 0) and vanishing nearest-neighbor interaction 
V = . This corresponds to the case of a uniform Hub- 
bard model which has been discussed in Refs. @, E3- 

In the upper and lower panels of Fig. [2] we plot respec- 
tively the electron density and the exact effective poten- 
tial for U = 6 and wo = 8. We consider in particular the 
case when the bulk is half-filled, i.e. N = Lb, and the 
cases of a bulk above and below half-filling, N = Lb =t 1. 

As shown in the lower panel, while the change in the 
potential on the bulk on going from N = Lb — 1 case to 
the TV = Lb is small and can be attributed to a 0(1 /TV) 
effect, there is a sizable [0(1)] jump on going from N = 
Lb to N = Lb + 1. For all other fillings different from 
N = Lb we find that the one particle addition jump is 
[0(1/N)\. 

The jump for N = Lb determines A xc which we define 
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FIG. 2: Panel (a) and (b) show respectively the charge den- 
sity and the KS potential for U = 6 at three different fillings, 
namely N = Lb, N = L^il. The inset of panel (a) presents 
a logarithmic plot of the density in the vacuum while the in- 
set of panel (b) presents a schematic comparison between the 
spectrum of the Hubbard model and of the effective KS sys- 
tem for N = Lb + 1. Other parameters: wo = 8t, v x = 0, 
Lb = 6 and Ly = 11. 
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tJ N+l 
J KS 
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V KS 
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with v^ s the TV-particle Kohn-Sham potential and 
ip^f (x), the corresponding z/-th eigenvector and eigen- 
value respectively. Notice that the shift of the potential 
is not perfectly constant in all the bulk region due to fi- 
nite size effects. Assuming that (p^^_\(x) is bound in the 
bulk region Eq. (j3j) correctly converges to the expected 
constant shift of the bulk in the thermodynamic limit. 
For finite systems we show below that with the present 
definition Eq. (pQ) is satisfied with surprisingly small fi- 
nite size corrections. Similar results are obtained if the 
slightly different definition of Ref. [6] is used [39]. 

Fig. [3] shows the U dependence of the exact charge gap 
of the N = Lb electron system defined as Ac = I N — A N 
where I N and A N indicate respectively the ionization 
energy and the electron affinity of the TV-particle sys- 
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E^ 1 with 



denoting the ground state energy of the TV-particle sys- 
tem obtained with the same Lanczos computation. We 
also show A xc + Aks where Aks is the exact KS gap, i.e. 
the gap in the spectrum of the effective non-interacting 
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TV-particle Kohn-Sham system. We see that indeed Eq.Q] 
is well fulfilled. As discussed below the Kohn-Sham gap 
should vanish in the thermodynamic limit for a Hubbard 
chain so its flniteness is a finite size effect. 

The charge density in the vacuum remains for all fill- 
ings much smaller than 1 and decays exponentially as 
shown by the logarithmic plot in the inset of Fig. [2] (up- 
per panel). The change in the density decay rate in the 
vacuum as the filling becomes larger than one (N > Lb), 
reflects a change in the ionization energy due to electronic 
correlations. Indeed as explained e.g. by Almbladh and 
von Barth HI, HH the density decay rate, ft, is related 
to the ionization energy. In particular in the lattice one 
can show that [Hj n = 2 cosh -1 (l N /2t). An accurate 
computation of the density profile in the vacuum region is 
what allow us to compute the absolute value of the Kohn- 
Sham potential in the bulk. More precisely the potential 
in the bulk is referred to the vacuum site furthest to the 
interface which is assumed to have zero Kohn-Sham po- 
tential. 

The inset of Fig. [3] shows schematically the behav- 
ior of the Kohn-Sham bands in a large Hubbard chain 
which can be solved exactly with periodic boundary 
conditions [40]. The charge is uniform and thus the Kohn- 
Sham potential is a constant which, without the vac- 
uum, remains undetermined. However we know that the 
chemical potential as a function of filling has a jump at 
half- filling equal to the Mott- Hubbard gap Amou- If we 
loosely consider the atoms of the Hubbard chain to have 
a large constant binding energy v x = — wo and to be 
immersed in a "vacuum" with zero binding energy we 
expect that the ionization energy will have a jump at 
half-filling due to the jump in the chemical potential. 
Due to DFT Koopman's theorem the center of the bands 
and Kohn-Sham potential will have a jump at half- filling 
such that A xc = Amou as shown schematically in the 
inset of Fig. [3l In Fig. [3] we also see that in spite of the 
bulk chain being short (Lb = 6), A xc approximately co- 
incides with A Mott fc> r the infinite system calculated by 
Bethe Ansatz showing that this picture @, EtJ * s indeed 
correct and finite size corrections to A xc are negligible. 

In Figure [2] we also note the appearance of a peak at 
the boundary between vacuum and bulk, on the vacuum 
side with width and height depend on the filling. Just 
as the leading (smallest) decay rate of the wave func- 
tion is determined by the first ionization energy, ioniza- 
tion from deeper states will determine subleading decays 
rates which are important at short distances [18|- Thus 
to better understand the origin of this peak it is useful 
to compare the photoemission spectrum of the Hubbard 
model and the Kohn-Sham spectrum. At large U and 
for N = Lb + 1 particles the removal spectra of both 
systems is very different as shown schematically in the 
inset of Fig. O In the Hubbard model only two states 
are available at low energy [42| while in the Kohn-Sham 
spectrum we have Lb + 1 states available. This large 
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FIG. 3: Exact charge gap Ac, Kohn-Sham gap, Aks and 
contribution of the xc-potential jump, A xc for a half- filled 
Hubbard chain with Lb — 6 sites. A Mott is the Mott gap 
for an infinite system calculated using Bethe Ansatz 40] . The 
inset shows Kohn-Sham band structure of a uniform Hubbard 
chain at half-filling (N — Lb) and with one added electron 
(JV = L B + 1). 



spectral difference would imply different subleading de- 
cay rates, with a tendency of the Kohn-Sham system to 
have a charge density larger than in the interacting sys- 
tem close to the boundary. This tendency is compensated 
by the appearance of the peak in the Kohn-Sham po- 
tential. Thus the anomalous transfer of spectral weight 
in the Hubbard model, which is the hallmark of strong 
electron correlation [42], reflects in the appearance of the 
barrier. 

Now we consider the transition between a Mott in- 
sulator an a ionic insulator. In order to simulate a bi- 
nary compound we consider the case in the presence of 
5, the Hubbard U which for simplicity is taken equal 
on all atoms and a nearest neighbor repulsion V . The 
system shows a transition form a ionic insulating regime 
to a Mott insulating regime when U ~ 25 + zV with 
z = 2 the coordination number [36|. In the atomic limit 
one finds that A^ ott = U — 26 in the Mott regime and 
j^iomc _ 2S -f 2zV — U in the ionic regime with both 
gaps coinciding at the transition. Notice that the latter 
is larger than the nearest neighbor charge transfer en- 
ergy corresponding to the excitation of a Frenkel exciton 
A ex = A^ omc — V and which becomes relevant bellow. 

Figure [H shows again that that Eq. (Q]) is well satisfied 
with negligibly finite size corrections [38]. (a) and (b) 
show respectively the results for V = and 5 = 2t and 
for V = O.bt and S = t. As one can easily check the total 
charge gap at U = for small t is the same in the two 
cases. However in the first case we have Ac ~ Aks in 
the ionic insulator and Ac — A xc in the Mott-insulating 
phase, while in the second case we have a finite contri- 



4 




2 4 6 8 




2 4 6 8 

U/t 



FIG. 4: Panel (a) and (b) show the different contributions to 
the gap, Aks and A xc and compare their sum to the exact 
charge gap calculated by Lanczos diagonalization, Ac- The 
parameters in the two panels are chosen to have the same total 
charge gap at U=0 in the small hopping limit. Parameters are 
in panel (a) S = 2£, V = 0, in panel (b) 8 = t, V = 0.5£. In 
both panels we set wo = —6t + U/2 and the potential of the 
site closer to the bulk- vacuum boundary has been chosen to 
correct boundary effects. 



bution of A xc to the gap in both regimes. Clearly the 
appearance of a finite A xc in the ionic regime is linked to 
the presence of the non local interaction V. This can be 
easily understand by considering the limit of weak tunnel- 
ing t << A ex . Using perturbation theory one easily finds 
that the amount of charge transferred from odd to even 
sites is 5p = 4t 2 / ' A 2 ex . For a uniform chain, by symmetry, 
the difference in the Kohn Sham potential between even 
and odd sites is equal to the Kohn Sham gap. Apply- 
ing the same perturbative argument to the Kohn- Sham 
system we arrive to the conclusion that to match the ex- 
act density = A ex therefore A xc = V. It is easy 
to check that these relations are valid in any dimension. 
They are in good agreement with the numerical results 
of Fig. [4] in the ionic regime. 

In general we expect that in strongly ionic insulators to 
a good approximation the Kohn-Sham gap matches the 
first Frenkel exciton and that A xc is given by its binding 
energy respect to the fundamental gap. While in ionic 
salts the Frenkel exciton is easily accessible experimen- 
tally the fundamental gap is difficult to measure and is 
often obtain by a theoretical fit to the observed optical 



spectra [43|. In any case matching of the Frenkel gap by 
puts a strong constraint on density functionals in 
strong ionic insulators. 

To conclude we have computed the exact exchange cor- 
relation potential of a correlated extended system includ- 
ing the (usually undetermined) absolute value respect to 
a vacuum level. This has allowed the first explicit nu- 
merical test of Eq. (pQ) in a model ionic/Mott insulator 
which dissipates any possible doubt on the validity of this 
equation or of the underling DFT-Koopmans theorem. 
For Mott insulators we have shown that the discontinu- 
ity of the exchange correlation potential is given by the 
Mott Hubbard gap which is of the order U for strong cor- 
relation. On the other hand in a strong ionic insulator 
the discontinuity is determined by the nearest neighbor 
repulsion V which provides a simple estimate of this elu- 
sive quantity. In addition we have shown that a surface 
correlation barrier appears in the effective potential of a 
correlated system when the removal spectrum of the sys- 
tem is very different from the removal spectrum of the 
Kohn-Sham system as is expected to occur in electron 
doped Mott insulators. 

This work was supported by the Italian Institute of 
Technology through the project NEWDFESCM. V.B. is 
indebted to L. Chiodo for discussions. 
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Exact Exchange correlation potential and proof of 
Koopmans theorem in the lattice 

To calculate the exact Kohn-Sham (KS) potential, i.e. 
the effective non-interacting potential which corresponds 
to the exact density, we adopt the following strategy: we 
first obtain the ground-state density by applying Lanczos 
diagonalization 4 to the "bulk+vacuum" chain, we then 
extract the KS potential by minimizing the difference 
between the KS density and the exact one for all values 
of the KS potential. 

The KS density is as usual expressed in terms of KS 
orbitals, cpi(x,a), as 



The orbitals tpi(x, a) are in turn defined through the well- 
known KS equations, 



(f 8 + v K s[x; p])ipi(x, a) = enpi(x, a) 



(1) 



where T s is defined by T s (pi(x,a) = —t(tpi(x + l,cr) + 
(fi(x — 1,<j) — 2ipi(x,a)), vks[x',p] is the effective KS 
potential and ei are the KS energies. To find the exact KS 
potential, we thus simply minimize relative mean square 
error on the density i.e. we calculate: 



mm 

v K s{x) 



\p(x,cr) - Pks(x,ct)\ 2 



\p(x,a) 



(2) 



where p denotes the exact density obtained by Lanczos 
diagonalization. After the minimization the relative er- 
ror on the density is smaller than 10 -5 i.e. |1 — pxs/p\ ^ 
10 -5 . As we now show, such a high accuracy is necessary 
to correctly describe the asymptotic decay of the density 
in the vacuum and to satisfy "Koopmans theorem" of 
DFT. This theorem, which identifies the highest occu- 
pied KS eigenvalue with the exact ionization energy, has 
been discussed and proved by several authors for stan- 
dard DFT in the "continuum" , here we extend the proof 
to lattice systems 

To prove "Koopmans theorem" in the lattice we fol- 
low the route underlined by Almbladh and von Barth 3 
and we introduce the quasiparticle amplitudes, f 3 (x) = 
(N — 1, s\c xa \N, 0), where |JV,0) and |JV, s) denote re- 
spectively the ground and the 8-th excited state of the 
full Af-electron Hamiltonian. 

Considering that in the vacuum sites the density is very 
small and decays exponentially with the distance from 
the surface for our parameter choice then, following the 




FIG. 1: Ionization energy computed from the exact ground 
state energies compared with the highest occupied Kohn- 
Sham eigenvalue for the system with Lb and Lb + 1 particles 
as a function of U/t for constant bulk potential wo = —St. 



same reasoning as in Ref. 3, one can show that many- 
body effects become asymptotically irrelevant. There- 
fore, deep in the vacuum, the quasiparticle amplitudes 
satisfy the Schrodinger like equation 

- t{f x+1 , a + /*_ 1>(T - 2/* CT ) - f^AE s = (3) 

where AE S indicates the difference between the ground 
state energy of the A^-particle system and energy s- 
th excited level of the TV — 1-particle system, AE S = 
Eq 1 — E^ -1 . Solving this recursive equation we obtain 
the following result for the decay of the ground-state 
quasi-particle amplitude, fo in the vacuum: 



fo(x) 



-(3x 



(4) 



with (3 = cosh -1 (4^), In being the ionization energy 
In = —AEq of the Af-particle system. For comparison 
in the continuum the decay rate is given by /o( r ) oc e~ Kr 
with k = \/2m e I where m e denotes the electron mass. 

Now considering that the density can be expressed 
in terms of the quasi-particle amplitudes as, p^ a = 
J2 S fs( xa )fs( xo ')i an d that all the quasi-particle ampli- 
tudes with s > decay exponentially faster than /o, we 
can assume that the exponential tail of the density far 
from the surface of a finite system will be governed by 
the ground state amplitude /q. The information on the 
ionization energy is therefore encoded in the decay of 
the density far from the bulk system's surface and it can 
be obtained in a DFT calculation. In particular, the 
relation between the ionization energy and the highest 
Kohn-Sham eigenvalue becomes clear if we consider the 
analogy between equation ([3|) and Kohn-Sham equations, 
Eqs. (pQ). We note indeed that, since the effective poten- 
tial can be chosen to vanish in the vacuum far from the 
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surface, in Kohn-Sham formalism the decay of the density 
is simply governed by the highest occupied Kohn-Sham 
level which thus concides with the ionization energy. 

The ionization energy as a function of U is shown in 
Fig. Q] • we see that Koopmans theorem of DFT is very 
well verified, the small error, within 2%, is due to the 
finite size of the vacuum chain. 



Charge gap, Kohn-Sham gap and exchange 
correlation potential discontinuity 

For sake of completeness we recall here the derivation 
of the relation between the charge-gap and the exchange- 
correlation potential shift starting from Koopmans the- 
orem, this also will allow us to identify the corrections 
due to finite size effects. 

By applying Koopmans theorem Ac can be rewritten 
as follows: 



JV+l 
-JV+l 



■ e N+l 



(5) 



where indicates the iV-th Kohn-Sham level of the M- 
particle system and the gap in the Kohn-Sham spectrum 



N 

AT+1 



of the iV-particle system is given by = e 

Using Kohn-Sham equations for the system with N 
and N + 1 particles the charge gap can be further parti- 
tioned as the sum of four terms, namely 



A c = A KS + A a 



(6) 



where A xc is the exchange-correlation gap introduced in 
the article, essentially equivalent to the one used in the 
literature (see e.g. Ref.[l|) 

A xc = \<P%tK x )\ 2 «s +1 (*) - <sW) (7) 

X 

while A^ and At are due to the relaxation of the iV+l-th 
Kohn-Sham orbital induced by the addition of a particle, 

X 

(8) 



A v = E (l^tl(^)l 2 - l^+i (*)| 2 ) <sW- (9) 



where v^ s (x) denotes the effective potential of the TV- 
particle system. 

For a translationally invariant system both A^ and 
At become of order 1/N and can be neglected in the 
thermodynamic limit. Since the difference Svks(x) = 
v ks^ 1 ( x ) — v ks( x ) becomes constant, in this limit, equa- 
tion [6] reduces to the well-know relation between the 
charge gap and the exchange-correlation potential dis- 
continuity shown e.g. in Ref. [2]. 

Interestingly we find that A^ and At tend to cancel in 
finite systems in all interaction regimes. As an example 




2 4 6 

U/t 



FIG. 2: Corrections to the charge gap, and At as func- 
tions of the ratio, U/t. The bulk potential is v x = 2(— l) x 
while the nearest-neighbor interaction parmeter V equals 
zero. Similar results are obtained both for V ^ and for 
v x - 0. 



in Figure [2] we plot the two contributions as a function of 
U/t for the ionic Hubbard model with v x = 2(— l) x and 
V = 0, as one can see the two terms are finite but with 
opposite signs. 

By applying perturbation theory it is not difficult to 
show that the cancellation stems directly from the defi- 
nition of exchange-correlation gap given in Eq. (0). The 
effective potential, v^ 1 , °f the + 1 particle system 



can be indeed related to the potential, v^ s , as follows 



..iV+l 



vg s + C + A*; 



(10) 



where C is a constant shift, C = 1/{Lb + 
Lv) ^2 X (^ks^ 1 — v ks)i which accounts for the disconti- 
nuity and remains finite in the infinite size limit and Av is 
a weak site dependent modulation such that Av = 
which vanishes in the infinite size limit. Starting from 
Eq. ([TO]) , by applying perturbation theory in Av we ob- 
tain the following expression for £^+1 



JV+l 
Z N+1 



£ N+1 



C + (<p% +1 (x)\Av\<f% +1 ) + 



E 

i/#iV+l 



\{<p» +1 \Av\<p») 



-N 
~"N+1 



(11) 



Similarly, by inserting Eq. (flQ|) in Eq. © we can recast 

^^nfC clS 



C+(^+l|A^|^+l) 



(12) 



JV+l 
-AT+1 



-AT+1 



Eventually, considering that A^ + At : 

A xc and expanding the wave- function \<p%X\) in Ec i- ® 
we arrive at the following result: 



E 



\(<p% +1 \Av\<p?) 



-JV+l 



riV 



(13) 



3 



From this equation we see that the sum A^ + A T is of 
second order in Av, since the change in the wave- function 
is of first order in A v. 

With a similar reasoning one can prove that with the 
definition of A xc given in Ref. [1] the finite-size errors to 



second order in Av have the same modulus but opposite 
sign. 
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